data = data3d('../data/head.npz', 'short')
step = float_param(0.002, 0.0001, 0.01, 'S')
threshold = float_param(-400, -800, 800, 'T')
cut = float_param(-1.8, -1.8, 1.8, 'C')
ambient = float_param(0.2, 0, 1, 'A')


# ray direction
D = normalize(E - S)

#light position
light_pos = float3(-1, -1, 1)

ray_length = length(E - S)
t = 0.0

while t < ray_length:
    P = S + t * D
    t += step
    
    f = cubic_query_3d_cut(data, P, D, cut) * 32768 - threshold
    if f < 0: continue
    
    Q = P
    P -= step * D
    
    # use bisection to refine the intersection between points P and Q
    while length(P - Q) > 0.000001:
        M = (P + Q) / 2
        f = linear_query_3d_cut(data, M, D, cut) * 32768 - threshold
        if f < 0:
            P = M
        else:
            Q = M
            
    # light direction
    L = normalize(light_pos - P)
    
    # normal vector
    N = -normalize(cubic_gradient_3d_cut(data, (P + Q) / 2, D, cut))
    
    # phong shading
    return phong(L, N, -D, float3(0.8, 0, 0), 0.2, 20, ambient)
    
return 0
